Modeling the role of the thalamus in resting-state functional connectivity: Nature or structure

The thalamus is a central brain structure that serves as a relay station for sensory inputs from the periphery to the cortex and regulates cortical arousal. Traditionally, it has been regarded as a passive relay that transmits information between brain regions. However, recent studies have suggested that the thalamus may also play a role in shaping functional connectivity (FC) in a task-based context. Based on this idea, we hypothesized that due to its centrality in the network and its involvement in cortical activation, the thalamus may also contribute to resting-state FC, a key neurological biomarker widely used to characterize brain function in health and disease. To investigate this hypothesis, we constructed ten in-silico brain network models based on neuroimaging data (MEG, MRI, and dwMRI), and simulated them including and excluding the thalamus, and raising the noise into thalamus to represent the afferences related to the reticular activating system (RAS) and the relay of peripheral sensory inputs. We simulated brain activity and compared the resulting FC to their empirical MEG counterparts to evaluate model’s performance. Results showed that a parceled version of the thalamus with higher noise, able to drive damped cortical oscillators, enhanced the match to empirical FC. However, with an already active self-oscillatory cortex, no impact on the dynamics was observed when introducing the thalamus. We also demonstrated that the enhanced performance was not related to the structural connectivity of the thalamus, but to its higher noisy inputs. Additionally, we highlighted the relevance of a balanced signal-to-noise ratio in thalamus to allow it to propagate its own dynamics. In conclusion, our study sheds light on the role of the thalamus in shaping brain dynamics and FC in resting-state and allowed us to discuss the general role of criticality in the brain at the mesoscale level.


Introduction
In humans, the thalamus is a nut size structure near the center of the brain that relays sensory inputs traveling to the cortex [1], fosters cortico-cortical communication through transthalamic pathways [2,3], and controls cortical arousal through the reticular activating system (RAS) [4]. To carry out these tasks, it contains three functionally distinct parts [5]: dorsal, ventral, and intralaminar. The dorsal part communicates bidirectionally with the cortex establishing two schemes of information exchange [6,7] (see Fig 1): first-order relay, in which the thalamus receives subcortical and sensory inputs (i.e., driving inputs), relays them to the cortex through excitatory thalamocortical cells and gets back modulatory feedback from layer 6 pyramidal neurons (i.e., modulatory inputs); and higher order relay, in which the thalamus receives inputs from layer 5 pyramidal cells of a cortical region and relays them to another location in the cortex, creating a transthalamic pathway for cortico-cortical connections [2,3]. In the ventral part, the reticular nucleus' inhibitory neurons establish connections both with each other and with neurons in the dorsal nuclei, to regulate and foster communication inside the thalamus [8,9]. The interactions between the dorsal and ventral parts of the thalamus allow for the generation of sustainable oscillations of neural activity, such as delta oscillations, sleep spindles, and slow waves, that may propagate to the cortex and influence its dynamics [10][11][12][13][14][15].
The intralaminar part is involved in the RAS [16] delivering cholinergic and monoaminergic neurotransmission diffusely to the cortex and controlling arousal [4].
Through these pathways, the thalamus connects to a widespread set of cortical regions [18-20], playing a role in different psychological processes such as sleep [10,21,22], pain [23,24], memory and learning [25], attention [26][27][28], motor and sensory processing [9,15,29], and consciousness [30][31][32][33][34]. This role of the thalamus has been classically depicted as a passive relay that transfers information between brain regions [25,35], however recent findings are challenging this view [20,25,[36][37][38]. Specifically, Schmitt et al. [26] showed that the mediodorsal thalamus was amplifying a FC pattern supporting the representation of specific task-related rules. All the thalamocortical mechanisms described above (i.e., the dorsal relay of peripheral sensory inputs, the transthalamic pathways for cortico-cortical communication, and the RAS neurotransmission that activates the cortex) may contribute to define the FC in the brain. Our study is aimed at understanding how.
FC is defined as a correlation between spatially distant neurophysiological signals [39] representing the functional integration of psychological processes in distributed brain networks [40]. Early neuroimaging studies were focused on revealing the activity patterns underlying cognitive processes during task execution, using resting-state as a control condition [41]. However, later findings showed that the brain in resting-state has a rich intrinsic activity [42,43] related to automatic and unconscious cognitive processing [44,45]. Since then, resting-state FC (rsFC) has been used to characterize brain function in health and disease [46][47][48] usually considering it as a static measure. More recently, this approach has been extended to capture the temporal richness of the activity patterns in resting-state through the concept of dynamical FC (dFC; [49,50] which has been suggested to reflect ongoing cognitive processing and that may be more informative of brain function and dysfunction than the static form [51][52][53][54]. Both metrics support the characterization of healthy aging, for which a general decrease in static FC, complemented by a slowing and less complex dFC has been shown and related to changes in cognitive performance [55][56][57][58]. Interestingly, some authors have proposed that changes in the thalamocortical network may contribute substantially to the disruptions in FC and cognition during aging [59,60]. Understanding the mechanisms that underlie and control (d)FC is an important research question that is still undisclosed, especially in aging. Here, we hypothesize that a similar thalamic mechanism that has been shown to be involved in defining FC during task execution [26] might also be active in resting-state.
Computational modeling allows for the generation of in-silico versions of real brains and personalized brain dynamics [61] employing brain network models (BNM). A BNM is based on: a structural connectivity (SC) network derived from diffusion-weighted MRI that captures how brain regions are wired together, and a set of neural mass models (NMM) that reproduce the electrophysiological dynamics of each brain region. A widely studied NMM is the Jansen-Rit (JR; [62]), a biologically-inspired model of a cortical column that implements excitatory and inhibitory subpopulations to produce oscillatory activity. This model shows a bifurcation over a parameter representing the strength of its inputs [63,64] that will be used in our work to reproduce different modes of the thalamocortical interaction. In a system, a bifurcation occurs when a change in the value of a parameter (i.e., bifurcation parameter) produces a qualitative change in the behavior of the system. For the JR model, the bifurcation separates two different states: a fixed point state, where the model behaves as a damped oscillator (prebifurcation), and a limit cycle state, where the model autonomously oscillates (postbifurcation). These two states turned out to be relevant to understand our findings.
To investigate the potential contribution of the thalamus to rsFC, we built ten in-silico BNMs based on healthy subjects' neuroimaging data (MEG, MRI, dwMRI) using JR NMMs. We simulated them using: 1) three SC versions (i.e., parceled thalamus, pTh; thalamus as a single node, Th; without thalamus, woTh) to explore the effect of both the parcellation and the mere presence of cortico-cortical transthalamic pathways, and 2) implementing a higher noisy input to the thalamus to reproduce its participation in RAS and the presence of peripheral sensory relays. We compared the simulated FC and dynamical FC (dFC) to their empirical MEG counterparts to evaluate performance. Additionally, we performed further simulations to explore under which conditions the thalamus contributes to the rsFC, including a control experiment using the cortico-cerebellar network instead of the thalamocortical one, and a set of parameter explorations over the intrinsic thalamic oscillatory behavior and the magnitude of the implemented noise. Our results showed that a limited set of driving nodes leading cortical activity was a plausible scenario in rsFC, where the thalamus would play a major role due to its nature: involved in the RAS system, and projecting sensory relays. These results contribute to the understanding of the basic principles of whole-brain function in health and disease, and to enrich the current picture of criticality behavior in the brain.

The thalamus impacts rsFC through its afferences
To explore the role of the thalamus in rsFC, we used two features of our in-silico BNMs: its structure, by simulating three different SC versions per subject (pTh, Th, woTh), and the NMMs noisy input by implementing higher than cortex thalamic noise (η th = [0.022, 2.2x10 -8 ], η cx = [2.2x10 -8 ]) to represent thalamic RAS system and peripheral first-order relays. We used the coupling parameter (g) to scale connectivity weights looking for the best match to empirical rsFC (i.e., the working point), as usual in whole-brain modeling [65][66][67][68]. Given that g acts as a bifurcation parameter, the simulated activity could be categorized into two regimes: prebifurcation, in which nodes operate as damped oscillators, and postbifurcation, in which nodes operate as autonomous oscillators. We simulated 60 seconds of brain activity per model and measured FC and dFC in the alpha band to compare to their empirical counterparts using Pearson's correlation (r PLV(α) ) and Kolmogorov-Smirnov distance (KSD), respectively. We will show statistical comparisons using the best values of those metrics per subject and bifurcation side.
Results on r PLV(α) showed a significant impact for both the structure [F(2, 18) = 191.77, In prebifurcation, implementing high thalamic noise raised significantly r PLV(α) values from close to zero to r PLV(α) �0.33 for Th [W = 0, Cohen's d = 5.09, p-corr = 0.005], and r PLV(α) �0.45 for pTh [W = 0, Cohen's d = 6.8, p-corr = 0.005] (see Fig 2). The resulting correlation values with high noise differed significantly between the implementations of the thalamic structure [F(2, 18)=186, Z 2 g ¼ 0:87, eps = 0.53, p-corr<0.0001], and pTh showed a global peak of r PLV(α) that unexpectedly overcame the values observed in postbifurcation in 8 out of 10 subjects (see S1 Fig) Fig 2). Interestingly, we noticed that high values of KSD(α) for woTh and Th were due to opposite underlying dFC distributions. woTh correlations were centered near r = 0 implying that FC matrices in time changed randomly, while Th correlations were centered near r = 1 implying that FC matrices in time were quasistatic (see S2 Fig).
In summary, the inclusion of the thalamus in the model had an impact just in prebifurcation range (in which nodes operate as damped oscillators), and only when implementing a high thalamic noise. For pTh, this condition overcame the performance of any other model. In postbifurcation range (in which nodes autonomously oscillate) the values for r PLV(α) were also high, however, the thalamus did not show a significant impact. The slight differences observed in the boxplots may be related to a higher number of autonomous oscillators trying to impose Simulations for two levels of thalamic noise (low noise in A, B, C; high noise in E, F, G), three implementations of thalamocortical SC (in colours), and exploring the parameter space for coupling factor (g) divided into prebifurcation (g < 7 shadowed regions in B and F) and postbifurcation (g > 7). B, F shows group averaged r PLV(α) and KSD(α) metrics. In the margins, boxplots show maximum values for those metrics per subject in prebifurcation (A, E) and postbifurcation (C, G). D shows the averaged bifurcation diagrams consisting of the maximum and minimum voltages per simulation for the cortex (thick line) and the thalamus (dashed line).
https://doi.org/10.1371/journal.pcbi.1011007.g002 their own dynamics that may make it more difficult to establish stable functional interactions between the nodes.

Structure is not the key: Comparing thalamocortical and cortico-cerebellar networks
We wondered whether the observed improvement of r PLV(α) with high thalamic noise in prebifurcation could be explained by the specific characteristics of the SC pattern of the thalamus. To test this, we performed a control experiment comparing the thalamocortical network to another with similar properties: the cortico-cerebellar network (see Table 1). We built three SC versions per subject: parceled cerebellum (pCer), single node cerebellum (Cer) and without cerebellum (woCer), all of them including the thalamus parceled. We simulated them implementing high noise into the cerebellum to compare the model performance to the thalamocortical network.
We observed that the general trend found with the thalamocortical networks persisted. woCer simulations showed close to zero r PLV(α) values in prebifurcation range with high noise, while Cer and pCer increased significantly their maximum correlations [F(2, 18) = 148.39, Z 2 g ¼ 0:895, eps = 0.75, p-corr<0.0001] up to similar values obtained with the thalamocortical network (see Fig 3). Moreover, pCer in prebifurcation also resulted in a global maximum in r PLV(α) compared to postbifurcation values. Interestingly, Cer changed the underlying bifurcation diagram of the model, moving it toward higher values of coupling. This was reflected by the peak of r PLV(α) in higher g values.
These results suggest that the specific thalamocortical SC pattern is not a major determinant of the thalamic contribution to rsFC.

Brain dynamics underlying each scenario
To understand the dynamics that underlie the observed values of r PLV(α) and KSD, we extracted a simulation sample per model condition with pTh (i.e., high/low thalamic noise, and pre-/ post-bifurcation; see Fig 4).
The dynamics in the prebifurcation range with high and low noise showed 1/f pink noise pattern. This is the result of a damped JR node processing a Gaussian noise as shown in previous literature [69,70]. High and low noise conditions in prebifurcation could be differentiated by their spectral powers and by the differences in FC matrices: with low noise, nodes are not powerful enough to interact, producing a functional disconnection that was captured by the FC and dFC matrices. In sharp contrast, in postbifurcation, nodes were self-oscillating in alpha frequency around 10Hz. Note the similarity between high and low noise to the thalamus in postbifurcation range, supporting the results reported in previous sections. Averaged network metrics for the thalamus, the cerebellum, and the global average of all regions. Similar relation to average was observed between regions in all metrics (degree, betweenness, and path length) except for node strength. For more details on the network analysis, see S1 Table. https://doi.org/10.1371/journal.pcbi.1011007.t001

Thalamic alpha propagates to the cortex within a balanced SNR
The best-performing thalamocortical model was the one in prebifurcation that integrated high thalamic noise and pTh. Looking at its underlying dynamics, we observed that the spectrum was showing a 1/f shape. As we are trying to reproduce MEG FC, in which a predominant alpha frequency is usually observed, we wondered whether a spectral change towards alpha would impact the model performance. In this section, we manipulated the oscillatory frequency of the thalamus making it surpass bifurcation and self-oscillate in alpha by varying its average input, p th . These simulations were performed with pTh structure and η th = 0.022. Simulations rising p th in prebifurcation (g<7), resulted in a transition of thalamic nodes from the noisy 1/f pink noise spectra to an alpha oscillation (see Fig 5, FFT peak), passing first through the slow and high amplitude limit cycle of the JR model [63,64] at p th = [0.11-0.13] (see the dark orange spot in Fig 5 SNR). In this transition along the bifurcation, r PLV(α) lowered down right after the high power and slow limit cycle (p th > 0.13). This could be due to a high SNR � 6 producing hypersynchronization (mean PLV�0.7; see Fig 5 mean PLV). However, we did not find this phenomenon with the slow limit cycle (PLV(α)mean �0.5) even though it showed a higher SNR. This might be explained by the alpha band filtering that leaves out of the analysis the potentially hypersynchronizing oscillations of the slow limit cycle.
Interestingly, the changes in thalamic activity indirectly increased the inter-regional inputs to cortical nodes, making some of them pass bifurcation at g�3 and g�6 (see the horizontal blue lines in Fig 5, FFT peak). These nodes produced a further rise in r PLV(α) (see the horizontal red line in Fig 5, r PLV(α) ) suggesting two important things for prebifurcation simulations: 1) that every node is a potential contributor of a general driving mechanism that we have located in the thalamus (through a high noise), and 2) that the number of drivers participating in that mechanism matters. In the range of p th = [0.13-0.3], where we observed alpha oscillations, r PLV(α) decreased. We wondered whether this effect could be related to the rise in SNR after setting the thalamus to oscillate in the alpha band. Therefore, we fixed p th = 0.15, and we varied the noisy input to the thalamus (η th ) to explore its impact on model performance. We found an optimal balance for SNR at η th = [0.05-0.15] in which the noisy inputs to the thalamus were enough to avoid hypersynchronization (see Fig 6 mean PLV) and low enough to maintain the intrinsic dynamics produced by the thalamus (i.e., alpha oscillations, see Fig 6 FFT peak). Further increases in noise would replace progressively the alpha oscillatory behavior by a 1/f spectra without reducing r PLV(α) .
From these observations, it could be thought that a general rise of noise in the model (i.e., to all nodes) would lead to better performance, however, in our modeling framework this is only true when the noise is implemented into a limited number of nodes. Independent noise into all cortical nodes is not linked to an enhancement of r PLV(α) (see S3 Fig).
From these parameter explorations, we extracted four additional models of interest to explore their underlying dynamics (see Fig 7). Two of them related to the exploration of p th : one for the hypersynchrony situation (p th = 0.15, η th = 0.022), another for the slow JR limit cycle (p th = 0.12, η th = 0.022); and another two regarding the exploration of SNR: one for the optimal SNR (p th = 0.15, η th = 0.09), and another for a higher noise than the optimal range (p th = 0.15, η th = 0.5). Fig 7 shows the underlying dynamics for each of the additional model scenarios simulated with g = 2.

Discussion
Understanding thalamocortical networks is crucial for unraveling the complex dynamics of the human brain. These networks are essential for transmitting sensory information from the periphery to the cortex and regulating cortical arousal, which are fundamental processes for perception, attention, and cognition. In this study, we aimed to gain a deeper understanding of the role of the thalamus in rsFC by utilizing computational brain models. Our experiments focused on testing two key features of thalamocortical networks: the presence of thalamic afferences related to the RAS and first-order sensory relays, and the structure of the thalamus. Interestingly, we found that only when we raised thalamic noisy inputs to represent the RAS and peripheral afferences activating a damped cortex, its presence affected the simulated dynamics. To validate our findings, we performed a control experiment using the cerebro-cerebellar network and showed that implementing high noise into the cerebellum could replicate the results observed in the thalamocortical network. Finally, we explored how the oscillatory Heatmaps showing different metrics from the same each simulation: the empiricalsimulated correlation of PLV in alpha band (r PLV(α) ), the simulated mean and std of the PLV(α) values, the frequency peak of the nodes' averaged spectra (FFT peak), the signal to noise ratio in thalamic nodes computed as the amplitude of simulated signals divided by the standard deviation of the Gaussian noise used for the thalamus (SNR(th)), and the bifurcation of cortical signals using the averaged maximum-minimum signals' voltage.
https://doi.org/10.1371/journal.pcbi.1011007.g005 behavior of thalamic nodes, specifically their frequency and SNR, could shape the emergent rsFC. Our study provides novel insights into the role of thalamocortical networks in shaping brain dynamics and highlights the relevance of balanced SNR activity for the propagation of alpha rhythms from the thalamus.
We expected that introducing the thalamus in our simulations would have generated a difference in model performance in postbifurcation, where the best correlations are usually found. However, this was not the case, as implementing the thalamus in our model only affected results in prebifurcation and when introducing a high thalamic noise to represent afferences from RAS system and peripheral sensory relays. In prebifurcation, nodes are operating as damped oscillators tending to relaxation at a fixed point. When we introduced high noise in the thalamus, its afferences to cortex rose cortical activation levels allowing for functional interactions. In this situation, the thalamus is driving cortical activation. Additionally, the best model performance was observed within this parametrization and including the thalamus with its dorsal nuclei divisions (pTh). More importantly, the postbifurcation range has been often associated with a state of generalized epileptic seizure [71,72] due to its highly synchronized intra-node oscillatory activity (see the spectra at postbifurcation in Fig 4B) that is not the empirical target of this study (i.e., MEG resting state). Taken together, these ideas led us to focus on the prebifurcation range.
A further control experiment demonstrated that this mechanism was not directly related to thalamic SC. This was ascertained by embedding the driver mechanism into another hub-like region (i.e., cerebellum) and obtaining equivalent results for the match to empirical rsFC. This would mean that in our model, where all neural masses are functionally equal, any node could theoretically play the role of the thalamus. In this line, during parameter explorations, some isolated cortical nodes passing bifurcation would contribute to enhance performance by being part of that driving mechanism. This would add up to the observation that parceled structures (i.e., thalamus and cerebellum) performed better than their single node versions, indicating that the number of driver nodes matters. In further research, it should be explored whether there is a computational optimum number of drivers to reproduce rsFC that could be subject or session-specific, and whether these differences may be linked to differences in brain and cognitive functioning. Taken together, our main conclusion is that a limited set of driving nodes is likely to underlie the dynamics of rsFC. We believe that those drivers might be linked both to the regions participating in the RAS system (Intralaminar Thalamus, Raphe Nuclei, and Locus Coeruleus) [4,16] and to the dorsal nuclei of the thalamus that are implicated in the relay of sensory information and have also been tightly linked to the generation of oscillatory behavior in the cortex in slow waves [73][74][75][76]. This would support the view of the thalamus as a driver and controller of cortical dynamics [6,26,[77][78][79].
Further explorations on the spectral characteristics of the drivers showed that the thalamus could propagate its own intrinsic alpha dynamics when a balanced SNR was achieved. This feature represents the interaction between thalamocortical pacemaker neurons [14] and peripheral sensory inputs reaching the system and provoking event-related desynchronization [80,81]. The model showed a good performance for reproducing empirical rsFC in that optimal range and with additional noise, generating a 1/f spectra. However, lower levels of noise with an alpha-oscillating thalamus reduced performance and led to a hypersynchronization situation transmitted from the thalamus to the cortex, generating an epileptic-like dynamic [82,83]. In line with these results, the thalamus has been proposed to be involved in the onset of temporal lobe epileptic seizures, transmitting more regular patterns of activity to the hippocampus [84,85].
From the computational perspective, previous work has shown that BNMs may reproduce better empirical rsFC when the models operate at the edge of bifurcation [67], the critical point. At this point, noisy excursions or the effect of nodes' interaction can lead the masses to behave in any of the two regimes separated by the critical point. This phenomenon, referred to as criticality [66], has been proposed to enhance the capacity of brain systems to convey information [86]. However, in our study, we showed an equivalent performance both at the edge of bifurcation and over the whole prebifurcation range. This contrast might be explained by the different cortical and thalamic parametrization implemented in our nodes. Our thalamic driving nodes feed the cortex, leading the dynamics. In the cited study [67], the nodes that randomly switch between states would represent the same driving mechanism as in our model. Interestingly though, our approach would suggest that criticality is not a necessary feature in the dynamics of a resting brain at the mesoscale level.
Many studies in the field have explicitly [87][88][89] or implicitly [66,[90][91][92] excluded subcortical regions from their BNMs. This could be due to the technical limitations of recording deep brain signals and/or to the complexity of reconstructing SC schemes for small, deep crossing-fibers regions. But, more importantly, it could be due to the unknown role that these regions may play in shaping simulated whole brain dynamics. Some studies have attempted to unravel these mechanisms, showing the importance of the cerebellum for brain dynamics [68], the relevance of cortico-subcortical interaction for shaping dynamical functional connectivity [93] and the relevance of neurotransmission [65,94]. Additionally, other studies are paving the way towards multiscale computational models in which subcortical areas are implemented with a further spatiotemporal level of detail [95][96][97][98]. These approaches could be an interesting path to extend our knowledge regarding the potential role of the thalamus (and other activating brain regions) in whole brain simulations.
In conclusion, our study provides novel insights into the role of thalamocortical networks in shaping brain dynamics. We demonstrate that a limited set of driving nodes leading cortical activation may better describe resting-state activity. The thalamus would be a relevant part in this mechanism due to its participation in the RAS system and through its peripheral sensory relays being delivered from its multiple dorsal nuclei. In this type of architecture, driving nodes might show a balanced SNR to avoid hypersynchronization in the network. Although it is still debated whether the thalamus has an active role in cognition [6,25,99], our study strongly suggests its active participation in driving cortical dynamics and shaping FC in resting-state. These findings may contribute to a better understanding of brain function and dysfunction, fostering the development of new therapeutic approaches targeting thalamocortical circuits.

Empirical dataset
MRI (T1 and DWI) scans and MEG recordings were acquired from 10 healthy participants in resting-state, with ages between 62 and 77 years old (mean 69, sd 4.17, 3 males, 7 females) from a dataset owned by the Centre of Cognitive and Computational Neuroscience, UCM, Madrid.
MRI-T1 scans were recorded in a General Electric 1.5 Tesla magnetic resonance scanner, using a high-resolution antenna and a homogenization PURE filter (fast spoiled gradient echo sequence, with parameters: repetition time/echo time/inversion time = 11.2/4.2/450 ms; flip angle = 12˚; slice thickness = 1 mm, 256×256 matrix, and field of view = 256 mm).
MEG recordings were acquired with an Elekta-Neuromag MEG system with 306 channels at 1000Hz sampling frequency and an online band-pass filtered between 0.1 and 330Hz. MEG protocol consisted of 5 min resting-state eyes closed.
All participants provided informed consent.

Functional connectivity
MEG recordings were preprocessed offline using the spatiotemporal signal space separation (tSSS) filtering algorithm [100], embedded in the Maxfilter Software v2.2 (correlation limit of 0.9 and correlation window of 10 seconds), to eliminate magnetic noise and compensate for head movements during the recording. Continuous MEG data were preprocessed using the Fieldtrip Toolbox [101], where an independent component-based algorithm was applied to remove the effects of ocular and cardiac signals from the data, together with external noise. Source reconstruction was performed using the software Brainstorm [102], anatomically informed by the MRI scans of each subject. We employed the minimum norm estimates method [103], with the constrained dipoles variant, by which the current dipoles are oriented normally to the cortical surface, to model the orientation of the macrocolumns of pyramidal neurons, perpendicular to the cortex [104] Source-space signals were then filtered in the alpha band (8)(9)(10)(11)(12) to calculate FC between the time series using the Phase Locking Value (PLV(α), [105]), and the resulting matrices were averaged into the AAL2 parcellation scheme [106]. We restricted the analysis to 1) cortical regions, to avoid the limitations of MEG recordings regarding deep brain signals [107]; and 2) the alpha band, for computational simplicity and being aware that it dominates MEG restingstate recordings. In addition, we computed dynamical functional connectivity matrices by extracting PLV(α) on consecutive intervals of 4 seconds of length with the sliding window approach and 50% of overlapping [108], and evaluating the correlation between these PLV(α) matrices.

Structural connectivity
Diffusion-weighted images were processed using DSI Studio (http://dsi-studio.labsolver.org). The quality of the images was checked before fiber tracking and corrected for motion artifacts, eddy currents, and phase distortions. Then, tensor metrics were calculated. To improve reproducibility, we used a deterministic fiber tracking algorithm with augmented tracking strategies [109][110][111]. The whole brain volume was used as seeding region. Both the anisotropy and angular thresholds were randomly selected (the latter, from 15 degrees to 90 degrees). The step size was randomly selected from 0.5 voxels to 1.5 voxels. A total of 5 million seeds were placed and tracks with lengths shorter than 15 or longer than 180 mm were discarded.
To explore the impact of including/excluding the thalamus in simulations, we performed a first experiment comparing three different structural connectivity (SC) versions of each subject brains': woTh, Th, and pTh. pTh consists of a brain network with 148 regions extracted from AAL3 atlas from which we kept the thalamic parcellation and, removed and merged the other areas to make it comparable to AAL2 scheme; Th consists of the 120 regions from AAL2; same for woTh in which we removed thalamic nuclei (118 regions). These three versions of the structural connectomes are represented in Fig 8 and included the cerebellum parceled into its nuclei. A list with all ROIs included can be found in S1 Table. Two connectivity matrices were calculated per SC version: counting the number of tracts connecting (i.e., passing through) each pair of brain regions, and the average length of those tracts.
As a control experiment, we applied the same process to the cerebellum as it is a brain region with similar network characteristics (see Table 1 and S1 Table), and it can also be modeled as a parceled structure and a single node. We extracted three SC versions per brain using AAL3 atlas: parceled cerebellum (pCer), cerebellum as a single node (Cer), and without cerebellum (woCer). The thalamus was modeled as parceled through all these versions, and therefore the resulting SCs were composed of 148, 122 and 120 regions, respectively.

Brain network model
SC matrices served as the skeleton for the BNMs implemented in TVB [112] where regional signals were simulated using JR NMMs [62]. This is a biologically inspired model of a cortical column capable of reproducing alpha oscillations through a system of second-order coupled differential equations: _ y 3 i ðtÞ ¼ AaS½y 1 i ðtÞ À y 2 i ðtÞ� À 2ay 3 i ðtÞ À a 2 y 0 i ðtÞ ð4Þ Where: The inter-regional communication introduces heterogeneity in terms of connection strength w ji , and conduction delays d ji (i.e., tract length / conduction speed) between nodes i and j, where: It represents the electrophysiological activity (in voltage) from three subpopulations of neurons: pyramidal neurons (y 0 ), excitatory interneurons (y 1 ), and inhibitory interneurons (y 2 ). These subpopulations are interconnected (Fig 9) and integrate external inputs from other cortical columns. The communication is implemented in terms of firing rate (Eqs 1 to 6) and a sigmoidal function (Eq 7) stands for the conversion from voltage to firing rate.
The input represents two main drivers of activity in the NMMs: inter-regional communication and intrinsic input. The former consists of the signal transmission between nodes through the SC of the brain in which weights are linearly scaled by a global coupling factor g, and tract lengths are divided by conduction speed to define d ji . Conduction speed was set to 15 m/s given the low impact shown in previous parameter space explorations done in this project (see S5 Fig). The latter is defined by a Gaussian noise with p mean and η std. Parameter values are described in Table 2.
The JR model shows two supercritical hopf bifurcations for the parameter p [64]. When JR NMMs are implemented in a connected network, the parameter g scales the inter-regional input to nodes, becoming a bifurcation parameter. We used the first bifurcation to separate two NMM's behaviors (Fig 10): damped oscillator in the prebifurcation range where nodes  In the center, a bifurcation diagram shows the minimum and maximum voltage for each value of the bifurcation parameter. At the critical point (dashed line), the bifurcation occurs and separates two states of the system: a) damped oscillator, whose activity tends to decay to a fixed point; and b) limit cycle oscillator, whose activity is a self-sustained oscillation. https://doi.org/10.1371/journal.pcbi.1011007.g010

PLOS COMPUTATIONAL BIOLOGY
The role of the thalamus in functional connectivity tend to a fixed point in voltage, and limit cycle in the postbifurcation range where nodes selfoscillate in the alpha frequency.

Simulations
In the first experiment, simulations were performed varying two parameters: the standard deviation of the input to the thalamus (i.e., noise, η th = [0.022, 2.2x10 -8 ]) representing the presence/absence of subcortical and peripheral inputs; and the thalamic structure (woTh, Th, pTh). The higher noise level was determined following previous research [12], while the lower was chosen to avoid flat signals in the fixed state of the JR model. Additionally, we explored the parameter space for coupling factor (g=[0-60]). These models were simulated with the parameter p set to 0.09. We performed three simulations of 60 seconds (removing the initial 4 seconds to avoid transients) per model and computed two metrics: the Pearson's correlation coefficient between the vectorized upper triangular matrices of the simulated and empirical FC (i.e., r PLV(α) ); and the KSD(α) between the distributions of correlations in the dFC matrices (empirical and simulated). The same configuration was used in the control experiment with the cerebellum. For further explorations, we simulated for 10 seconds (omitting the initial 2 seconds to avoid transients) different ranges of the mean intrinsic input to the thalamus p th ) and its standard deviation (η th ). Besides the r PLV(α) and KSD(α) metrics, we show 1) bifurcation diagrams capturing the averaged maximum and minimum signal's voltage per simulated ROI at each point in the parameter space. In the case of exploring 2 parameters at the same time (e.g., g and p th ), bifurcations are presented in heatmaps conveying information about the difference between the maximum and minimum signal voltage; 2) signal-to-noise ratio (SNR) in the thalamus that is computed by dividing the amplitude of simulated signals by the standard deviation of the Gaussian noise used for the thalamus; 3) relative power between cortex and thalamus calculated by dividing the averaged area of cortical spectra by the averaged area of thalamic spectra.

Statistics
We averaged the results of the 3 sets of simulations (i.e., repetitions) and performed statistical analysis for the group of 10 subjects. For the first experiment, we considered the maximum r PLV(α) and minimum KSD(α) per subject, thalamic SC version and scenario. The effects of thalamic SC version and noise levels were evaluated using four two-way repeated measures ANOVA: two comparing maximum r PLV(α) in prebifurcation and postbifurcation; and another two comparing minimum KSD; after checking for the statistical assumptions of normality (Shapiro's test) and sphericity (Mauchly's test). Pairwise comparisons for thalamic structure and thalamic noise were evaluated using Wilcoxon test, correcting for multiple comparisons using FDR Benjamini-Hochberg method. The same procedure was applied in the control experiment with the cerebro-cerebellar network to maximum r PLV(α) comparisons. in both prebifurcation (g = [3,7]) and postbifurcation 2 (g = [9,40]). The three thalamocortical SC (i.e., woTh, Th, pTh) versions were simulated.  Table. Network analysis of the regions included in the BNMs. Degree, the number of neighbors of a region, and node strength, the average number of streamlines connecting a region to others were normalized over their respective maxima. Betweenness captures the number of shortest paths in a network that passes through a node. Path length stands for the average of the shortest paths for a node. Metrics were calculated with Networkx package in Python 3.9. The thalamus and the cerebellum are considered here in the parceled version. Note that Cingulate_Ant in AAL3 is divided in 3 parts and it was merged to match AAL2 scheme.